!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!> \file
!! Simple ODE program to test the numerical integrators.
!!
!! \n
!!
!! This example has two purposes: set a common structure for the time
!! integrators and provide a specific example to test them.
!!
!! The building block for any time integrator is the type
!! <tt>t_time_state</tt>, which is used to represent the state of the
!! solution and is described in details in \c mod_time_integrators.
!<

module ode_test_problem

 use mod_messages, only: &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   wp

 use mod_time_integrators, only: &
   c_ode

 implicit none

 public :: &
   t_ode, ode_exact

 private

 type, extends(c_ode) :: t_ode
  real(wp) :: f1(4)
  real(wp), allocatable :: f2(:)
 contains
  procedure, pass(tnd) :: rhs => ode_rhs
  procedure, pass(z)   :: add => ode_add
  procedure, pass(y)   :: stimes => ode_stimes
  procedure, pass(x)   :: solve  => ode_solve
 end type t_ode

 integer, parameter :: &
   nn     =  1
 real(wp), parameter :: &
   lambda = -3.0_wp, &
   tau    =  7.3_wp, &
   alpha  =  0.1_wp

contains

 subroutine ode_rhs(tnd,t,uuu)
  real(wp),     intent(in)    :: t   !< time level
  class(c_ode), intent(in)    :: uuu !< present state
  class(t_ode), intent(inout) :: tnd !< tendency

   select type(uuu); type is(t_ode)

   tnd%f1(1) = lambda*uuu%f1(1)
   tnd%f1(2) = lambda*uuu%f1(2)**2
   tnd%f1(3) = lambda*sign(1.0_wp,uuu%f1(3))*sqrt(abs(uuu%f1(3)))
   tnd%f1(4) = alpha * (t**nn)

   tnd%f2(1) = -tau*uuu%f2(2)
   tnd%f2(2) =  tau*uuu%f2(1)

   end select
 end subroutine ode_rhs

 pure subroutine ode_add(z,x,y)
  class(c_ode), intent(in)    :: x, y
  class(t_ode), intent(inout) :: z

   select type(x); type is(t_ode); select type(y); type is(t_ode)

   z%f1 = x%f1 + y%f1
   z%f2 = x%f2 + y%f2

   end select; end select
 end subroutine ode_add

 pure subroutine ode_stimes(y,r,x)
  real(wp),     intent(in)    :: r
  class(c_ode), intent(in)    :: x
  class(t_ode), intent(inout) :: y

   select type(x); type is(t_ode)

   y%f1 = r*x%f1
   y%f2 = r*x%f2

   end select
 end subroutine ode_stimes

 pure subroutine ode_solve(t,x,sigma,b,xl)
  real(wp),     intent(in) :: t, sigma
  class(c_ode), intent(in) :: b, xl
  class(t_ode), intent(inout) :: x

  real(wp) :: v1, v2

   select type(b); type is(t_ode); select type(xl); type is(t_ode)

   x%f1(1) = 1.0_wp/(1.0_wp-lambda*sigma) * b%f1(1)
   x%f1(2) = 1.0_wp/(1.0_wp-lambda*sigma*xl%f1(2)) * b%f1(2)
   v1 = sqrt(abs(xl%f1(3)))
   v2 = lambda*sigma/2.0_wp
   x%f1(3) = v1/(v1-v2)*( b%f1(3) + v2*sign(1.0_wp,xl%f1(3))*v1 )
   x%f1(4) = b%f1(4) + sigma*alpha * (t**nn)

   x%f2(1) = b%f2(1) - sigma*tau*b%f2(2)
   x%f2(2) = sigma*tau*b%f2(1) + b%f2(2)
   x%f2 = x%f2/(1.0_wp+(sigma*tau)**2)

   end select; end select
 end subroutine ode_solve

 pure subroutine ode_exact(t,y)
  real(wp),     intent(in) :: t
  class(t_ode), intent(inout) :: y

   y%f1(1) = 3.0_wp*exp(lambda*t)
   y%f1(2) = 2.0_wp/(1.0_wp-lambda*2.0_wp*t)
   y%f1(3) = max(lambda*t/2.0_wp+1.0_wp,0.0_wp)**2
   y%f1(4) = 0.01_wp + alpha/real(nn+1,wp) * (t**(nn+1))

   y%f2(1) = cos(tau*t)
   y%f2(2) = sin(tau*t)

 end subroutine ode_exact

end module ode_test_problem


program ode_test

 use mod_messages, only: &
   mod_messages_constructor, &
   mod_messages_destructor,  &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_constructor, &
   mod_kinds_destructor,  &
   wp

 use mod_time_integrators, only: &
   mod_time_integrators_constructor, &
   mod_time_integrators_destructor,  &
   ee_init,      ee_step,      ee_clean,      &
   hm_init,      hm_step,      hm_clean,      &
   ssprk54_init, ssprk54_step, ssprk54_clean, &
   thetam_init,  thetam_step,  thetam_clean,  &
   bdf2_init,    bdf2_step,    bdf2_clean,    &
   bdf3_init,    bdf3_step,    bdf3_clean,    &
   bdf2ex_init,  bdf2ex_step,  bdf2ex_clean


 use ode_test_problem, only: &
   t_ode, ode_exact

 integer :: n, n0
 real(wp), parameter :: dt = 0.05_wp/2.0_wp
 real(wp) :: t
 type(t_ode) :: uuu0, uuun, work(8), uuuex

 ! Note: the following commands are useful to check the results in
 ! octave
 !  load ode-test.out
 !  figure;plot(ode_test(:,1),ode_test(:,8:13),'k');hold on;plot(ode_test(:,1),ode_test(:,2:7),'o-');figure;plot(ode_test(:,1),ode_test(:,14:end))
 !  sum(abs(ode_test(:,14:end)),1)*4/size(ode_test,1)

  call mod_messages_constructor()
  call mod_kinds_constructor()


  call mod_time_integrators_constructor()

  allocate(uuu0%f2(2),uuun%f2(2),uuuex%f2(2))
  call ode_exact(0.0_wp,uuu0)

  ! allocate also the work array
  do n=1,size(work)
    allocate(work(n)%f2(2))
  enddo
  open(22,file="ode-test.out", status='replace',action='write', &
       form='formatted',position='rewind')

  n0 = 0 ! should be changed for multistep methods
  call ee_init(dt)
  call hm_init(dt)
  call thetam_init(dt,0.5_wp)
  call ssprk54_init(dt)
  call bdf2_init(dt,0.0_wp,uuu0,work(3)); n0=1
  !call bdf3_init(dt,0.0_wp,uuu0,work(1:3),work(4:5)); n0=2
  !call bdf2ex_init(dt,0.0_wp,uuu0,work(1:2),work(3:4)); n0=1
  do n=n0,nint(4.0_wp/dt)-1
    t = real(n,wp)*dt

    !call ee_step(uuun,t,uuu0,work)
    !call hm_step(uuun,t,uuu0,work)
    !call thetam_step(uuun,t,uuu0,work)
    !call ssprk54_step(uuun,t,uuu0,work)
    call bdf2_step(uuun,t,uuu0,work(1:2),work(3))
    !call bdf3_step(uuun,t,uuu0,work(1:3),work(4:5))
    !call bdf2ex_step(uuun,t,uuu0,work(1:2),work(3:4))
    uuu0 = uuun

    call ode_exact(t+dt,uuuex)
    write(22,'(19e14.6)') &
      t+dt,               &
      uuu0%f1,  uuu0%f2,  &
      uuuex%f1, uuuex%f2, &
      uuuex%f1 - uuu0%f1, uuuex%f2 - uuu0%f2
  enddo
  call ee_clean()
  call hm_clean()
  call thetam_clean()
  call bdf2_clean()
  call bdf3_clean()
  call bdf2ex_clean()
  call ssprk54_clean()

  close(22)

  call mod_time_integrators_destructor()


  call mod_kinds_destructor()
  call mod_messages_destructor()

end program ode_test
